#实习七

#-------------practice1-------------------

#对ADdata采用log2进行标准化，并利用ANOVA检验分别求所有蛋白在3组之间的P值
rm(list = ls())
load("AData1.RData")
#将0替换成NA
match_merge_data[match_merge_data == 0] <- NA
Adata <- match_merge_data
group <- as.data.frame(match_merge_data[,101])
colnames(group) <- "group"
#log2标准化
Adata <- log2(Adata[,1:100])
Adata <- cbind(group,Adata)
pos_ad <- group == "ad"
pos_asym <- group == "asym"
pos_ctl <- group == "ctl"
pvalue <- c(1:100)
for(i in 2:101){
  #统计每一种类型的非NA值个数
    ad <- sum(!is.na(Adata[pos_ad,i]))
    ctl <- sum(!is.na(Adata[pos_ctl,i]))
    asym <- sum(!is.na(Adata[pos_asym,i]))
    #每个组中非NA值均大于三个时利用oneway.test计算p值
    if(ad >= 3 & ctl >= 3 & asym >= 3){
      p <- oneway.test(Adata[,i]~group, Adata)
      pvalue[i-1] <- p$p.value
    }
    #否则将p值赋NA
    else{
      pvalue[i -1] = NA
    }
}
#以数据框形式并且合并蛋白名字一起输出p值列表
protein <- as.data.frame(colnames(Adata))
protein <- as.data.frame(protein[-1,])
pvalue <- as.data.frame(cbind(protein,pvalue))
colnames(pvalue) <- c("protein","P")
save(pvalue, file = "pvalue.RData")
pvalue

#------------practice2----------------------
#提取class6_volcano.RData数据中p<0.05的蛋白进行GO富集分析

#安装包clusterProfiler，org.Hs.eg.db
#参考链接1：http://www.bio-info-trainee.com/7479.html
#参考链接2：https://zhuanlan.zhihu.com/p/423710750
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install(version = "3.14")

options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor/")
options("repos" = c(CRAN="http://mirrors.cloud.tencent.com/CRAN/")) 
options(download.file.method = 'libcurl')
options(url.method='libcurl')

#利用BiocManager安装所需的R包
BiocManager::install("clusterProfiler",ask = F,update = F)
BiocManager::install("org.Hs.eg.db",ask = F,update = F)


#清空环境
rm(list = ls())

#加载所需包
library(clusterProfiler)
library(org.Hs.eg.db)
library(ggplot2)
#加载数据
load("class6_volcano.RData")
#数据预处理，筛选p值小于0.05的蛋白
P <- prostat$P
index <- P < 0.05
pro_data <- prostat[index,]


allID <- pro_data$ID[!is.na(pro_data$ID)]
enall <- enrichGO(allID,"org.Hs.eg.db",keyType = "SYMBOL",ont = "ALL",pvalueCutoff = 0.05,
                  pAdjustMethod = "BH",minGSSize = 10,maxGSSize = 500)

#条形图展示结果
barplot(enall, showCategory = 20,drop = T)

#气泡图展示结果
dotplot(enall, showCategory = 30)




#--------------------practice3-------------------
#对题2富集分析结果进行气泡图展示

rm(list = ls())
library(clusterProfiler)
library(org.Hs.eg.db)
library(ggplot2)
#加载数据
load("class6_volcano.RData")
#数据预处理，筛选p值小于0.05的蛋白
P <- prostat$P
index <- P < 0.05
pro_data <- prostat[index,]

allID <- pro_data$ID[!is.na(pro_data$ID)]
enall <- enrichGO(allID,"org.Hs.eg.db",keyType = "SYMBOL",ont = "ALL",pvalueCutoff = 0.05,
                  pAdjustMethod = "BH",minGSSize = 10,maxGSSize = 500)


#气泡图展示结果
dotplot(enall, showCategory = 30)
    